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The statistical properties of the dissipation process constrain the analysis of large scale numerical 
simulations of three dimensional incompressible magnetohydrodynamic (MHD) turbulence, such 
as those of Biskamp and Miiller [Phys. Plasmas 7, 4889 (2000)]. The structure functions of the 
turbulent flow are expected to display statistical self-similarity, but the relatively low Reynolds 
numbers attainable by direct numerical simulation, combined with the finite size of the system, make 
this difficult to measure directly. However, it is known that extended self-similarity, which constrains 
the ratio of scaling exponents of structure functions of different orders, is well satisfied. This implies 
the extension of physical scaling arguments beyond the inertial range into the dissipation range. 
The present work focuses on the scaling properties of the dissipation process itself. This provides 
an important consistency check in that we find that the ratio of dissipation structure function 
exponents is that predicted by the She and Leveque [Phys. Rev. Lett 72, 336 (1994)] theory 
proposed by Biskamp and Miiller. This supplies further evidence that the cascade mechanism 
in three dimensional MIfD turbulence is non-linear random eddy scrambling, with the level of 
intermittency determined by dissipation through the formation of current sheets. 
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I. INTRODUCTION 

This paper investigates the previously under explored topic of scaling in the local rate of dissipation in magne- 
tohydrodynamic (MHD) flows. Turbulent fluids and plasmas display three properties that motivate development 
of statistical theories (i) disorganisation, in the sense that structures arise on all scales; (ii) unpredictability of 
detailed behaviour, in the sense of inability to predict a signal's future behaviour from knowledge of its past, implying 
links with deterministic chaos; and (iii) reproducibility of statistical measures, combined with the presence of statisti- 
cal self-similarity. Much progress has been made by the heuristic treatment of scaling laws derived from energy cascade 
arguments, following Kolmogorov and Richardson, see for example Ref.jJ]. The basic idea is that energy-carrying 
structures (eddies) are injected on large scales, non-hnear eddy interaction causes energy to cascade to smaller scales 
in a self-similar manner, and energy is finally dissipated by viscosity on small scales. A quasi-stationary state evolves 
where the rate of viscous dissipation matches the rate of energy injection. Scaling exponents Cp characterise the 
' resulting statistical self-similarity found in structure functions Sf: 

Sf = ((v(x + 1, <) . l/l - v(x, t) ■ l/lf) - l^^ (1) 

Here v is the fluid velocity, x is a position vector, 1 is a differencing vector, and the average is an ensemble average. 
The statistical self-similarity represented by the power-law in I is only valid within the inertial range Id I lo', 
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here Iq is the characteristic macroscale, and Id is the dissipation scale at which the cascade terminates. The set 
of scaling exponents Cp in Eq.(^) is expected to be universal since it characterises the generic cascade process. It 
is worth noting here that universality can only be expected in the isotropic case. When anisotropics are present, 
deviation from the isotropic case can be expected, and this will relate to the strength of the anisotropy. In MHD 
turbulence, anisotropy can be introduced in the form of an imposed magnetic field. The effect of this on the scaling 
exponents is investigated in Ref.2]. This reference also investigates anisotropy in terms of that introduced by the 
local magnetic field even when an applied field in absent. This stems from the Goldreich and Sridhar objection to the 
assumption of local isotropy in the MHD case ^] . In Ref. ^ structure functions are calculated with the differencing 
length perpendicular and parallel to the local magnetic field. The perpendicular structure functions were found to 
exhibit stronger intermittency than the parallel structure functions. Exponents calculated from the perpendicular 
structure functions were found to coincide with those calculated from the isotropic structure functions. Essentially 
dimensional arguments, in which the relevant physical parameters are identified heuristically, have been formulated 
to provide basic fluid scaling information. These arguments linearly relate Cp to p, for example the Kolmogorov 1941 
phenomenology 0,0] predicts (p = p/3. As such, basic fluid scaling can be characterised by one number a such that 

Sf ~ (2) 

To exploit these concepts, let us write the equations of incompressible MHD in Elsasser symmetric form 0: 

dtz^ = -zT.Vz=^ - V + + {v/2 + r//2) V^z^ + {v/2 - ry/2) V^z^F (3) 



V.z± = (4) 

Here the Elsasser field variables are z* = v ± B (/iop) ^ , where p is the scalar pressure, v is kinematic viscosity, 
rj is magnetic diffusivity and p is fluid density. The symmetry of Eq.© suggests that statistical treatment of z^ 
may be more fundamental than separate treatments of v and B. In light of this, longitudinal structure functions are 
constructed in terms of Elsasser field variables hereafter: 

Sf ^^ = (I (z(±) (x + 1, i) • \/l - z(±) (x, t) ■ 1//) \P) ^ f'p^' (5) 

As mentioned above, heuristic arguments that make predictions about basic fluid scaling only linearly relate Cp to p. 
In reality depends nonlinearly on p due to the intermittent spatial distribution of eddy activity. Basic fluid scaling 
can be modified to take this into account by the application of an intermittency correction. A commonly applied class 
of intermittency correction describes statistical self-similarity in the local rate of dissipation ei by means of scaling 
exponents Tp: 




For a review of the fractal nature of the local rate of dissipation for hydrodynamics, see for example Ref. 



As we shall see, the intermittent nature of the system is captured by the nonlinear dependence of Tp on p in 
Eq.(|SJ|. This nonlinearity can conveniently be expressed multiplicatively in relation to the basic linear fluid scaling 
of Eq. Q . Specifically we may write 

Sf ^ {efr^ (7) 

Where g is a constant whose value is inferred from model assumptions such as those of Kolmogorov (K41) and 
Iroshnikov-Kraichnan (IK) 8, 9]. This is Kolmogorov's refined similarity hypothesis j^. The scaling exponents Cp 
in Eq.lO are inferred by Eq.JTJ to be Cp = Tpg + pa- That is, the intermittency in the velocity field structure 
functions is achieved via reasoning concerning local rate of dissipation. One model that uses this hypothesis, and has 
proven successful in predicting the scaling exponents for hydrodynamic turbulence, is that from the 1994 paper of 
She and Leveque (SL) [l^. Here physical assumptions are made regarding the scaling of the local rate of dissipation. 
Specifically: the hierarchical nature of Eq.lO above, as expressed in Eq.(6) of Ref.^^l; the rate of dissipation by 
the most intensely dissipating structures is related only to the eddy turnover time as determined by the basic fluid 
scaling, as in Eq.(4) of Ref.jia]; and the space filling nature of the most intensely dissipating structures can be 
described by one parameter (their Hausdorff dimension) . These three assumptions can be combined to formulate a 
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second order difference equation for the scaling exponents Tp that has one non-trivial solution, as in Eq.(9) of Ref.Tcl. 
This solution can be formulated in terms of the two parameters: the co-dimension of the most intensely dissipating 
structures, C = D — dff, where D is the embedding dimension and is the Hausdorff dimension; and the basic fluid 
scaling number expressed by a in Eq.Q above. Following Ref. , we may write 

Tp^-{l-a)p + C-Cil-il-a)/C)P (8) 

This two parameter formulation follows that previously noted by Dubrulle , whose parameters A and P correspond 
to our (1 — a) and 1 — A/C respectively. The refined similarity hypothesis, as expressed in Eq.jTJ), is then invoked to 
obtain the following expression for the structure function scaling exponents Cp- 

= pa - (1 - a)pg + C-C{1-{1- a)/Cfs (9) 

Previously Elsasser field structure functions have been identified with an SL model of the type Eq.©, see, for example 
Refs. 12, 13]. In the present paper the refined similarity hypothesis for MHD is tested by applying a modified form of 
Eq.(jZ|), see Ea. H12|l . to the simulation data of Biskamp and Miiller. This provides an important consistency check for 
previous studies. Equation ©, which probes the multifractal properties of the local rate of dissipation, but does not 
rely on the refined similarity hypothesis, can also be tested directly against the simulation results, as we discuss below. 



Direct numerical simulations must resolve the dissipation scale Id so that energy does not accumulate at large 
wavenumbers, artificially stunting the cascade. Most of the numerical resolution is therefore used on the dissipation 
range, whereas it is only on scales much larger than that dissipative effects are negligible, and scaling laws of 
the type discussed arise. Thus high Reynolds number simulations with an extensive inertial range are currently 
unavailable. However, the principle of extended self-similarity (ESS) 0| can be used to extend the inertial range 
scaling laws into the range of length scales that is significantly affected by dissipation but still larger than Id- Instead 
of considering the scaling of individual structure functions, the principle of ESS involves investigating the scaling of 
one order structure function against another, on the assumption that 

- {sf^f'^'-' (10) 

Here it can be seen that any set of structure functions will satisfy this relation providing 

Sf ^ G(/)^^ (11) 

where G(l) can be any function of I which is independent of p. Here we use the notation of S. C. Emily et al. that 
these authors used to describe the general properties of generalised extended self-similarity, as expressed in Eq.(8) 
of Ref.jl^ - though generalised ESS is not discussed in the present paper. When the objective of ESS is to extend 
scaling into the dissipation range, G{1) can be rewritten as l"^ where G'{1) is introduced to accommodate the non 
constant flux of energy through length scales in the dissipation range. As such, G'{1) asymptotically approaches one 
as the length scale increases from the dissipative to the inertial range. 



The She-Leveque model as it has appeared so far, is only valid in the inertial range. Let us now discuss how 
this model can be interpreted in the framework of ESS. This problem has been tackled for hydrodynamic turbulence 
by Dubrulle for example. In that paper the explicit inclusion of I in the refined similarity hypothesis [Eq.© 
with g = a = 1/3 for hydrodynamic turbulence] is replaced by a generalised length scale, which is cast in terms of 
the third order structure function as expressed in Eq.(12) of Ref.[lJ]. This problem was addressed similarly by Benzi 
et al. where the scaling relation 

Sf - (ef) (^y-)"" (12) 

is explicitly formulated in Ref . Il^ . The appropriate relation between <^p and Tp is now — Tpg +pa(^i/a- Using this 
relation combined with Ea. Hll|) and Ea. p2|l it can be seen that (ef) must also have the form (ef) ~ G{IY^ . This 
implies ESS exists also in the local rate of dissipation, such that 

(ef) ^ {e1)^^^'^^) (13) 

It can then be seen that if a She-Leveque model of the general type Eq.Q is used to explain scaling exponents 
obtained via ESS, as expressed in Ea. HlU|) . then two consistency checks are appropriate. First Kolmogorov's refined 
similarity hypothesis should be satisfied in the form Ea. (|12|) . and second ESS should exist in the local rate of 
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FIG. 1: Extended self-similarity for the Elsasser field variable z'"^^ (order five against order three), compare Eg. lfOL for decaying 
MHD turbulence where structure functions are normalised by the total energy before time averaging. This normalisation reveals 
the same underlying scaling for points from different simulation times, as shown. After Biskamp and Miiller fl^ . 



dissipation as in Ea. H13|) . 

The present paper performs these checks for consistency for the simulation of Biskamp and Miiller T^]. Here 
the scaling exponents C^^' [see Eq.©] were investigated via direct numerical simulation of the three dimensional (3D) 
incompressible MHD equations, with a spatial grid of 512^ points The simulation is of decaying turbulence 

with initially equal magnetic and kinetic energy densities and u = r]. A fuller discussion of the numerical procedure 
is present in the next section. Since the turbulence decays with time, structure functions are normalised by the 
total energy in the simulation (kinetic plus magnetic) before time averaging takes place. Biskamp and Miiller [T^ 
extracted the ratios of scaling exponents Cp/Cs by ESS and directly determined i^p ~ 1. These exponents were found 
to match a variant of the She-Leveque 1994 model Eq.@ inserting Kolmogorov basic fluid scaling {g = a = 1/3) 
with the most intensely dissipating structures being sheet-like. (C = 1). Early investigations of this type assumed 
Iroshnikov-Kraichnan fluid scaling where the most intensely dissipating structures are sheet-like |1^ (see Ref. 
[isl l for a derivation of Cp^^ for this case). Sheet-like intensely dissipating structures are thought to exist in MHD 
turbulence because of the propensity of current sheets to form. We refer to Fig. 5 of Ref.^^l foi' isosurfaces of current 

density squared, and to Fig|21for isosurfaces constructed from the shear in the z'"*"^ field . Bothfi gures show 

the existence of two dimensional coherent structures; Fig[5]is more directly related to the analyses presented in the 
present paper. Basic Kolmogorov fluid scaling for Alfvenic fluctuations has been verified for low Mach number (~ 0.1) 
compressible 13, 20] and incompressible 12, 17] 3D MHD turbulence by power spectrum analysis, and by checking for 
scaling in the third order structure function such that C3 = 1. Extended self-similarity has also been utilised to extract 
ratios of scaling exponents related to an inverse cascade in laboratory plasmas [23]. In other work, a generalised 
version of this SL model has been applied to compressible flows where C is allowed to vary as a fitting parameter 
[T3L l2^. and in the case of Ref. 22] this dimension is interpreted as a function of the sonic Mach number. Figure 
^shows an example of the normalisation and ESS procedure for z^^-* structure functions from the data analysed here. 

The philosophy behind our investigation can now be summarised as follows. Given a simulation, the set of 
structure functions Sf can be calculated. These are expected to display statistical self-similarity as expressed in 
Eq.(|SI, where the scaling exponents give insight into the physics of the cascade process. The relatively low 
Reynolds numbers attainable by direct numerical simulation, combined with the finite size of the system, make this 
statistical self-similarity difficult to measure directly. However, it is found that extended self-similarity of the type 
expressed in Ea. (|10|) is well satisfied, allowing the ratio of scaling exponents Cp/Cs to be directly measured. There 
is a range of a priori views concerning these ratios, reflecting physical model assumptions. The ratios of scaling 
exponents recovered from ESS analysis of the 3D MHD simulation data are compared with these models, and the 
best fit is identified. Our investigation thus assists in validating the physical assumptions made in formulating the 
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FIG. 2; (Color online) Isosurfaces of sheet-like (2D) coherent structures of the squared gradient of the z'^-* Elsasser field variable 
from the 3D MHD turbulence simulation ofBiskamp and Miiller. 



currently favoured model, namely Eq.(P|with g = a = 1/3 and C = 1 giving = p/9 + 1 - (l/3)P/3. In particular, 
we confirm the existence of a specific type of extended self-similarity in the local rate of dissipation, with exponents 
given by Eq.® with a = 1/3 and C — 1 giving Tp — — 2p/3 + 1 — (1/3)^. We also show that Kolmogorov's refined 
similarity hypothesis, in the form Ea. (|12|l . is satisfied. 



II. NUMERICAL PROCEDURES 

The data analysed here stems from a direct numerical simulation of decaying isotropic turbulence (see Ref.p^ 
for additional details). The equations of incompressible MHD are written in terms of the vorticity, w = V x v, 
in order to eliminate the pressure variable. These are solved by a pseudospectral scheme (see, for example, 
Ref . [23^ ) . Time integration is carried out by applying the trapezoidal leapfrog method The aliasing error 

associated with this approach is reduced to discretisation error level by spherical truncation of the Fourier grid |26j | . 

The simulation domain comprises a periodic box in Fourier space with 512'^ points. Initially the fields have 
random phases and amplitudes ~ exp(— fc^/(2A:Q)) with ko — 4. The ratio of total kinetic and magnetic energy of the 

fluctuations is set to one. Cross helicity, which is proportional to Jy ^(z+)^ — (z^)^^ dV, is absent throughout the 
duration of the simulation. The mag netic hehcity, = ^JyA-BdV, where B = V x A, is set to 0.7 H^^^. Here 
~ E^^ /ko where E^^ is the energy in the magnetic field. The diffusivities u = r] = 4 x 10^* imply a magnetic 
Prandtl number Pr„i = v/rj oi one. 

The run was performed over ten eddy turnover times, defined as the time required to reach maximum dissi- 
pation when starting from smooth initial fields. Structure functions and moments of dissipation are calculated in 
time intervals of 0.5 between t — 3.5 and t — 9.5, during which the turbulence undergoes self-similar decay. 
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III. RESULTS 

In the present paper, the gradient squared measure {diz'f^'^)'^ is used as a proxy Q for the local rate of dissipation 
{diBj - djBiY r]/2 + {d,Vj + djV.y v/2. This proxy has recently been employed to study turbulence in the thermal ion 
velocity of the solar wind as measured by the ACE spacecraft "27], giving results consistent with those presented below. 
This is particularly interesting insofar as the solar wind remains one of the few accessible system of quasistationary 
fully developed MHD turbulence @ although we note that MHD intermittency studies have also been performed 
on the reversed field pinch experiment RFX [2^. Figure [3 shows isosurface plots of the gradient squared measure 
for the simulation of Biskamp and Miillcr. Two dimensional coherent structures dominate the image, suggesting the 
dimension parameter entering an SL model should equal two, as in the model employed by Biskamp and Miiller. 
Following Eq.®, statistical self-similarity in the dissipation measure is expressed as 




This proxy, which involves a one dimensional i nteg ration rather than the full 3D integration of Eq.©, facilitates 
comparison with related experimental MHD [2lll27| and hydrodynamic studies, and also offers benefits in 

computation time. 

The SL model adopted by Biskamp and Miiller predicts 

rf) =-2p/3+l-(l/3f (15) 

This is simply Eq.JHJ with a = 1/3 and C = 1. Gradients are calculated from the data of Biskamp and Miiller 
using a high order finite difference scheme, and the integral is performed by the trapezium method. Normalisation 
by the spatial average of viscous plus Ohmic rates of dissipation allows time averaging to be performed. Figure El 
shows an example of the ESS and normalisation procedure for xf'^'' order p = 5 against order p = 3. Statistical 
self-similarity is recovered with roll-off from power law scaling as I aproaches the system size. This roll-off behaviour 
at large I may be due to the finite size of the system, since a more extensive part of the simulation domain is 
encompassed by the spatial average (the integral over dV) as / increases in Ea. (|14|l . In FigO points identified with 
this roll-off are removed, and ratio of scaling exponents (Tp/ra) is calculated from the remaining points by linear 
regression. These ratios are shown in Fig0| No significent difference between the scaling recovered from z(+) and 
7,^~^ can be seen. This should be expected since no theoretical distinction needs to be drawn between z'^-* and z^"' 
for the vanishing values of cross helicity present in this simulation. The solid line in Fig^l shows the ratio predicted 
by Ea. (|15|l . in contrast to the dashed line which shows the ratio predicted by the SL theory for hydrodynamic 
turbulence [To|. Caution must be taken when calculating high order moments, since these are strongly influenced 
by the tails of their associated distributions. This can easily lead to bad counting statistics. The order p is only 
taken up to p = 6.5 for the dissipation measure (instead of p = 8 as for the Elsasser field structure functions |12j |) 
because of the extremely intermittent nature of the signal; large values affect the average [the angular brackets in 
Ea. (|14|l ] more as the order p increases. This effect is evaluated using a similar methodology to that in Ref. 30]. If a 
worst case scenario is imagined, where the average is strongly affected by one point in the signal, one would expect 
1/51 members of the spatial average in Ea. (|14|l to be directly affected by this point, where 51 is the grid spacing. 
We can then define an event as incorporating 1/61 members of the spatial average. It is found that ~ 5 percent 
of the average is defined by only ~ 10 events for order p = 6.5. This situation is of course worse for higher values of p. 

Plots were constructed in order to test Ea. H12|l . This involves taking the product of structure functions of 
the field variables and the dissipative quantities, in contrast to FigsD] and El Figures O and El show these plots for 
n = 1.5 and n = 2 respectively. The low order measure in FigElshows a relation that is nearly linear, with a gradient 
close to the ideal value of one, see Ea. (|12() . This is encouraging considering the deviation expected at the smallest 
and largest scales due to finite size effects. However, unlike the case in FigEl there may be some curvature across the 
range of the plot. The higher order measure in FigEldeviates from a linear scaling relation. We note that taking this 
test to high order involves the product of two quantities that have challenging counting statistics, plotted against a 
high order structure function. The deviation of the gradient seen in FigjSlfrom the ideal value of one is perhaps not 
surprising, because the constraints described above become stronger at high order. 
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FIG. 3: Extended self-similarity in the Elsasser field variable z'+^ gradient squared proxy for the local rate of dissipation (order 
five against order three), compare Eg. 1131 with the gradient squared proxy from Eg. 11411 replacing e^. Normalisation by the 
space averaged local rate of viscous and Ohmic dissipation allows time averaging in spite of the decay process. Deviation from 
power law scaling at large I is probably a finite size effect. The solid line is the best fit in the linear region. 
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FIG. 4: Ratio of scaling exponents (order p over order three) obtained via extended self-similarity from the Elsasser field 
gradient squared proxy for the local rate of dissipation. Errors in these measurements lie within the marker symbols. Solid line 
shows ratios predicted by a She-Leveque theory based on Kolmogorov fluid scaling and sheet-like most intensely dissipating 
structures, Eg. 11511 . The dashed line shows ratios predicted by hydrodynamic She-Leveque fTol| . 



IV. CONCLUSIONS 



Extended self-similarity is recovered in the gradient squared proxy for the local rate of dissipation of the Elsasser 
field variables z^^'^^ computed by Biskamp and Miiller. We believe this is the first time this has been shovi^n for MHD 
fiows. This result supports the application to Elsasser field scaling exponents (^p^^ of turbulence theories that require 
statistical self-similarity in the local rate of dissipation, even when (^p^^ are extracted from relatively low Reynolds 
number flows via ESS. Furthermore the ratio of exponents recovered is that predicted by the SL theory proposed 
by Biskamp and Miiller This supplies further evidence that the cascade mechanism in three dimensional MHD 
turbulence is non-linear random eddy scrambling, with the level of intermittency determined by dissipation through 
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FIG. 5: Plot to test Kolmogorov's refined similarity hypothesis as applied to extended self-similarity, Ea. p2|l . This involves 
taking the product of the field variable and dissipative structure functions, in contrast to Figs0and|3| Agreement with the 
hypothesis would give a straight line with unit gradient. Normalisation was performed as in FigsQandElto allow time averaging 
despite the decay process. 




FIG. 6: High order test of Kolmogorov's refined similarity hypothesis as applied to extended self-similarity. Eg. 1121 1. Normali- 
sation is performed as in Figl^ 



the formation of two dimensional coherent structures. However, Kolmogorov's ESS modified refined similarity 
hypothesis remains to be verified at high order. 
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